H3 Grid and Food Security#
This notebook showcases how the proposed H3 grids (level 6) relate to sub-national areas with accute food insecurity phases (IPC).
Here, we compare the latest data avialble for Kenya from FEWS NET and FAO/IPC, and propose an approach to merge the scores into the standardized grid.
Show code cell source
import sys, os
from os.path import join, expanduser, basename
import geopandas as gpd
sys.path.append(join("X:", 'Repos', 'pacific-observatory', 'scripts', 'python'))
import folium as flm
from grid_utils import *
from matplotlib.colors import ListedColormap
import ee
ee.Initialize()
from geemap.conversion import *
from geemap.foliumap import Map
from geemap import colormaps
DATA_DIR = join("X:", 'data', 'CRW')
target_path = os.path.join(DATA_DIR, 'Shapefile', 'FEWS_Admin_LZ_v3_simplified.shp')
adm0_path = os.path.join(DATA_DIR, 'Shapefile', 'FEWS_Admin_LZ_adm0.shp')
aoi = gpd.read_file(target_path)
aoi = aoi.to_crs('EPSG:4326')
countries = ['Kenya']
# # countries = ['Somalia', 'South Sudan', 'Haiti', 'Kenya', 'Afghanistan', 'Democratic Republic of the Congo']
aoi = aoi.loc[aoi.ADMIN0.isin(countries)].copy()
adm0 = gpd.read_file(adm0_path)
adm0 = adm0.loc[adm0.ADM0_NAME.isin(countries)].copy()
1. Latest IPC data from FEWS NET, Februrary 2023#
The FEWS NET data is provided at the intersection of sub-national areas and livelihood zones.
Show code cell source
fews_latest = join(DATA_DIR, 'notebook', 'true', '2023', 'may5', 'FEWS February 2023 Update TrueBoundaries_05-05-23.csv')
fews_df = pd.read_csv(fews_latest)
fews_df = fews_df.loc[fews_df.year_month=="2023_02"].copy()
fews_df = fews_df.loc[fews_df.country==countries[0]].copy()
aoi_df = aoi.merge(fews_df[['admin_code', 'fews_ipc_adjusted']], on="admin_code", how="left")
green = (220/255, 240/255, 220/255)
yellow = (250/255, 230/255, 30/255)
orange = (230/255, 120/255, 0/255)
red = (200/255, 0/255, 0/255)
darkred = (100/255, 0/255, 0/255)
cmap = ListedColormap([green, yellow, orange, red, darkred])
m = aoi_df.explore(
column='fews_ipc_adjusted',
tooltip= ["admin_name", "admin_code", "ADMIN0", "fews_ipc_adjusted"],
name="FEWS Admin LZ",
cmap = cmap,
categorical = True,
categories = [1,2,3,4,5],
legend_kwds=dict(caption='Integrated Phase Classification'),
legend=True
)
flm.LayerControl().add_to(m)
m
Make this Notebook Trusted to load map: File -> Trust Notebook
2. Map IPC data to H3 grid (level 6)#
We merge the IPC phases to the H3 grid using a spatial overlay, assigning phases based on largest overlap.
Show code cell source
adm0 = adm0.explode(index_parts=False)
h3_df = gdf_to_hex_df(adm0, 6, overfill=False)
h3_df.drop_duplicates('hex_id', inplace=True)
h3_gdf = h3list_to_gdf(h3_df)
h3_gdf_join = gpd.overlay(h3_gdf, aoi_df[['admin_code', 'admin_name', 'fews_ipc_adjusted', 'geometry']], how='intersection')
h3_gdf_join = h3_gdf_join.to_crs('EPSG:3857')
h3_gdf_join.loc[:, 'area'] = h3_gdf_join.area
h3_gdf_join.sort_values(by='area', inplace=True, ascending=False)
h3_gdf_join.drop_duplicates(subset='hex_id', keep='first', inplace=True)
h3_gdf_merge = h3_gdf.merge(h3_gdf_join[['hex_id', 'admin_code', 'admin_name', 'fews_ipc_adjusted']], on="hex_id", validate='1:1')
m = aoi_df.explore(
column='fews_ipc_adjusted',
tooltip= ["admin_name", "admin_code", "ADMIN0", "fews_ipc_adjusted"],
name="FEWS Admin LZ",
cmap = cmap,
categorical = True,
categories = [1,2,3,4,5],
legend_kwds=dict(caption='Integrated Phase Classification'),
legend=True
)
h3_gdf_merge.explore(
m = m,
column='fews_ipc_adjusted',
tooltip= ["admin_name", "admin_code", "fews_ipc_adjusted"],
name="H3 Grid Level 6",
cmap = cmap,
categorical = True,
categories = [1,2,3,4,5],
legend=False
)
flm.LayerControl().add_to(m)
m
Input polygon is too small to be filled in with hexagon of resolution: 6
Attempting to fit a hex of resolution: 7...
A hex of resolution: 7 fits the input polygon
Finding it's parent in resolution: 6
Make this Notebook Trusted to load map: File -> Trust Notebook
3. Compare FEWS NET and IPC data sources#
The following map shows differences in the assesments provided by FEWS NET and IPC for the same time period (Feb 2023).
Also note the difference in units, IPC provides data at the district level.
Show code cell source
gdf = gpd.read_file(join(expanduser("~"), 'Repos', 'food-security', 'ipc.geojson'))
gdf = gdf.loc[gdf.country=="KE"].copy()
m = aoi_df.explore(
column='fews_ipc_adjusted',
tooltip= ["admin_name", "admin_code", "ADMIN0", "fews_ipc_adjusted"],
name="FEWS NET",
cmap = cmap,
categorical = True,
categories = [1,2,3,4,5],
legend_kwds=dict(caption='Integrated Phase Classification'),
legend=True
)
gdf.explore(
m = m,
column='overall_phase',
tooltip= ["title", "id", "overall_phase", "country", "phase3_worse_percentage"],
name="IPC",
cmap = cmap,
categorical = True,
categories = [1,2,3,4,5],
legend=False
)
flm.LayerControl().add_to(m)
m
Make this Notebook Trusted to load map: File -> Trust Notebook
4. Overlay NDVI with H3 Grid#
In this map, we see how the grid relates to an NDVI raster from MODIS (250 mts.).
We plot both the raw NDVI raster, as well as a version masked with a global crop mask (GFSAD).
Show code cell source
# aoi_ee = gdf_to_ee(h3_gdf_merge.loc[h3_gdf_merge.admin_code==3076].copy())
aoi_ee = gdf_to_ee(h3_gdf_merge.loc[h3_gdf_merge.admin_code.isin([3089, 3077, 3076, 3082, 3062, 3074, 3086, 3070, 3091, 3128])].copy())
terra = ee.ImageCollection('MODIS/006/MOD13Q1').select('NDVI')
aqua = ee.ImageCollection('MODIS/006/MYD13Q1').select('NDVI')
modis = terra.merge(aqua)
crop_mask = ee.Image("USGS/GFSAD1000_V1").gt(0) #.clip(aoi_ee)
ken_ee = ee.FeatureCollection("FAO/GAUL_SIMPLIFIED_500m/2015/level0").filter(ee.Filter.eq('ADM0_NAME', 'Kenya'))
one = ee.Image(modis.first()).clip(ken_ee)
def applyMask(image, mask):
return image.updateMask(mask)
scale_factor = 0.0001
def applyScaleFactor(image, scale_factor):
return image.multiply(scale_factor)
one_scaled = applyScaleFactor(one, scale_factor)
one_masked = applyMask(one, crop_mask)
one_masked_scaled = applyScaleFactor(one_masked, scale_factor)
pal = colormaps.get_palette('BrBG', 10) # Greens
crop_vis = {
'min': 0,
'max': 1,
'palette': ['white', 'green'],
}
ndvi_vis = {
'min': 0.2,
'max': 0.6,
'palette': pal
}
centroid = aoi_ee.geometry().centroid().getInfo()['coordinates']
m = Map(center=(centroid[1], centroid[0]), zoom=9, add_google_map=False)
# m.setOptions(mapTypeId='HYBRID')
m.add_basemap('HYBRID')
# m.addLayer(crop_mask, crop_vis, 'Mask')
m.addLayer(one_scaled, ndvi_vis, 'NDVI')
m.addLayer(one_masked_scaled, ndvi_vis, 'NDVI Masked')
m.addLayer(aoi_ee, {}, 'H3 Grid')
m